LiMM-PCA with the Sensory Data
Manon Martin
July 19, 2020
1 Data set up
table(design\(Judges, design\)Candies)
| Â | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| 01 | 3 | 3 | 3 | 3 | 3 |
| 02 | 3 | 3 | 3 | 3 | 3 |
| 03 | 3 | 3 | 3 | 3 | 3 |
| 04 | 3 | 3 | 3 | 3 | 3 |
| 05 | 3 | 3 | 3 | 3 | 3 |
| 06 | 3 | 3 | 3 | 3 | 3 |
| 07 | 3 | 3 | 3 | 3 | 3 |
| 08 | 3 | 3 | 3 | 3 | 3 |
| 09 | 3 | 3 | 3 | 3 | 3 |
| 10 | 3 | 3 | 3 | 3 | 3 |
| 11 | 3 | 3 | 3 | 3 | 3 |
raw_outcomes <- outcomes
df <- data.frame(raw_outcomes, design)
library(doBy)
# Mean attributes values by judges
table_summary_mean <- summaryBy(. ~ Judges , data = df,
FUN = mean )
pander(table_summary_mean)| Judges | Transp.mean | Acid.mean | Sweet.mean | Raspb.mean | Sugar.mean |
|---|---|---|---|---|---|
| 01 | 9.29 | 4.69 | 7.83 | 4.43 | 5.75 |
| 02 | 8.36 | 3.69 | 4.2 | 4.49 | 5.73 |
| 03 | 9.95 | 4.78 | 6.99 | 7.07 | 5.7 |
| 04 | 9 | 3.22 | 4.28 | 3.73 | 6.1 |
| 05 | 5.96 | 7.59 | 9.27 | 7.92 | 4.72 |
| 06 | 8.67 | 5.55 | 7.02 | 5.73 | 4.81 |
| 07 | 7.64 | 10.36 | 7.68 | 5.99 | 5.19 |
| 08 | 8.26 | 4.17 | 6.96 | 4.2 | 5.67 |
| 09 | 8.29 | 5.79 | 7.69 | 6.07 | 6.04 |
| 10 | 8.76 | 4.92 | 3.98 | 5.11 | 5.98 |
| 11 | 8.35 | 3.39 | 7.77 | 6 | 5.31 |
| Bites.mean | Hard.mean | Elastic.mean | Sticky.mean |
|---|---|---|---|
| 10.08 | 9.21 | 8.54 | 8.87 |
| 6.48 | 7.14 | 7.93 | 7.33 |
| 9.12 | 8.83 | 7.53 | 7.55 |
| 9.66 | 9.57 | 9.32 | 7.26 |
| 8.5 | 8.96 | 9.02 | 7.08 |
| 9.13 | 9.46 | 8.79 | 8.89 |
| 9.39 | 9.26 | 9.45 | 9.49 |
| 8.88 | 8.36 | 7.75 | 8.79 |
| 9.7 | 9.35 | 9.1 | 9.41 |
| 9.5 | 8.68 | 9.68 | 8.29 |
| 9.58 | 9.53 | 8.58 | 7.79 |
2 remove outliers ??
exclude <- FALSE
if(exclude){
index <- which(!rownames(raw_outcomes) %in% c("0211", "0713"))
raw_outcomes <- raw_outcomes[index,]
design <- design[index,]
}2.1 Graph
# pdf(file.path(fig_path, "SDA_outcomes.pdf"), width = 8, height = 5,
# pointsize=19)
par(mar=c(4,4,2,5), xpd = TRUE)
col <- gg_color_hue(n=5)[design$Candies]
outcomes <- as.matrix(outcomes)
plot( outcomes[1,], type="l", xaxt="n",
ylim=range(outcomes), col=col[1],
xlab="Attribute", ylab = "Rating", main = "Sensory Data outcomes")
for (i in 2:dim(outcomes)[1]){
lines(outcomes[i,], col=col[i])
}
axis(side=1, at = 1:9, labels = colnames(outcomes))
legend("topright", legend=levels(design$Candies),
col=rainbow(n=5), lwd=1, title="Candies",
inset=c(-0.24,0), box.col = "white")
2.2 Hotelling T^2 for outlier detection
model = pca(raw_outcomes, scale = FALSE,
info = 'Simple PCA model', lim.type = "jm")
# model$ncomp.selected
ncomp <- 4
# plotVariance(model)
# plotResiduals(model, show.labels = TRUE, ncomp = ncomp,
# main = "Squared residual distance vs Hotelling T2 distance",
# norm = TRUE)
# plotResiduals(model, show.labels = TRUE, ncomp = ncomp,
# main = "Squared residual distance vs Hotelling T2 distance")
Qlim <- model$Qlim
T2lim <- model$T2lim
rownames(Qlim)[2] <- rownames(T2lim)[2] <- "Out_limit"
# In case of PCA the critical limits are just shown
# on residual plot as lines and can be used for detection
# of extreme objects (solid line) and outliers (dashed line).
plot_hotelling <- function(){
xlim <- range(model$calres$T2[,ncomp])
xlim[1] <- xlim[1]*0.9
xlim[2] <- xlim[2]*1.1
ylim <- range(model$calres$Q[,ncomp])
ylim[1] <- ylim[1]*0.9
ylim[2] <- ylim[2]*1.1
plot(model$calres$T2[,ncomp], model$calres$Q[,ncomp],
main = "Diagnostic plot for score \nand residual outliers", xlab = "Hotelling T2 distance",
ylab ="Squared residual distance", pch = 16, xlim = xlim,
ylim = ylim, cex.lab=1.25, cex.main=1.3)
abline(h=Qlim["Out_limit",ncomp], v=T2lim["Out_limit",ncomp], lty =3)
legend("topright", legend = "Outlier limit", lty = 3, cex=1.2)
index1 <- which(model$calres$T2[,ncomp]>=T2lim["Out_limit",ncomp])
index2 <- which(model$calres$Q[,ncomp]>=Qlim["Out_limit",ncomp])
index_ho <- unique(c(index1,index2))
text(x = model$calres$T2[index_ho,ncomp], y=model$calres$Q[index_ho,ncomp],
labels = names(model$calres$T2[index_ho,ncomp]), pos = c(1,2,3,4))
}
plot_hotelling()
# pdf(file.path(fig_path,"SDA_hotelling.pdf"), height = 6, width = 6,
# pointsize = 12)
plot_hotelling()
3 Distribution of the raw outcome variables
par(mfrow=c(1,3))
for (i in 1:ncol(raw_outcomes)){
# qqPlot(raw_outcomes[,i])
d = density(raw_outcomes[,i])
# hist(raw_outcomes[,i])
plot(d)
}
4 Dimension reduction by PCA
n = nrow(raw_outcomes)
###################################################
# Dimension reduction by PCA
###################################################
# ===== PCA ===== #
res_pca <- MBXUCL::SVDforPCA(x = raw_outcomes)
pander("Cumulated variance")Cumulated variance
| PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 | PC9 |
|---|---|---|---|---|---|---|---|---|
| 74.69 | 84.08 | 90.55 | 93.44 | 95.42 | 97.3 | 98.58 | 99.67 | 100 |
4.1 PCA plots
df <- data.frame(PC = as.character(1:length(res_pca$var)),
var = res_pca$var)
screePlot <- ggplot(df, aes(y=0,yend=var,x=PC,
xend=PC))+ geom_segment()+
labs(title= "Scree \nplot",
x = "PC", y="% var")+theme_classic()+
theme(axis.title.x=element_text(size=12),
axis.title.y=element_text(size=12))
index <- 27
pch <- c("1"=15, "2"=2, "3"=19, "4"=4, "5"=8)
scoresPlot <- DrawScores(res_pca, axes=c(1,2),drawNames = FALSE,
main = "Scores plot",
pch = Candies, size=2,
legend_shape_manual = pch) +
coord_cartesian(xlim = c(-25, 25), ylim = c(-15, 15)) +
annotate("text", y = (res_pca$scores[index,2] +
1*c(1)),
x = res_pca$scores[index,1],
label = rownames(res_pca$scores[,1:2])[index])## Scale for 'shape' is already present. Adding another scale for 'shape', which
## will replace the existing scale.loadPlot <- ScatterPlot(res_pca$loadings[,1], res_pca$loadings[,2],
points_labs = rownames(res_pca$loadings), cex.lab=4) +
geom_vline(xintercept = 0, lwd=0.1) +
geom_hline(yintercept = 0, lwd=0.1) +
xlab(label=paste0("PC1 (", round(res_pca$var[1],2),"%)")) +
ylab(paste0("PC2 (", round(res_pca$var[2],2),"%)")) +
coord_cartesian(xlim = c(-0.55, 0.55), ylim = c(-0.8, 0.8)) +
ggtitle(label = "Loadings plot")4.2 Apply first step of LiMMPCA (PCA)
5 Tests with univariate ANOVA
mat <- cbind(subject=1:dim(design)[1], design, outcomes)
mat$subject <- as.factor(mat$subject)
########################
# for one response
########################
# res_aov <- aov(PC1 ~ Candies*Judges + Error(subject/Judges), data=mat)
# sumar <- summary(res_aov)
# sumar <- sumar[[1]][[1]]
# Df Sum Sq Mean Sq F value Pr(>F)
# Candies 4 31470.6 7867.7 780.1762 < 2e-16 ***
# Judges 10 224.9 22.5 2.2304 0.02089 *
# Candies:Judges 40 707.7 17.7 1.7545 0.01158 *
# Residuals 110 1109.3 10.1
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# MeanSq <- sumar$`Mean Sq`
# names(MeanSq) <- gsub(" ", "", rownames(sumar))
# MeanSq["Candies"]
# MeanSq["Judges"]
# MeanSq["Candies:Judges"]
# MeanSq["Residuals"]
#
# MeanSq["Candies"]/MeanSq["Candies:Judges"]
# MeanSq["Judges"]/MeanSq["Residuals"]
# MeanSq["Candies:Judges"]/MeanSq["Residuals"]
########################
# for all the responses
########################
sumar_list <- vector(mode="list", length=nPC)
pvalues <- matrix(data = NA, nrow=nPC, ncol=3,
dimnames = list(NULL, c( "Candies","Judges","Candies:Judges")))
for (i in 1:nPC){
form <- paste(paste0("PC",i), "~ Candies*Judges + Error(subject/Judges)")
res_aov <- suppressWarnings(aov(as.formula(form), data=mat))
sumar_list[[i]] <- summary(res_aov)[[1]][[1]]
cat("PC ", i)
print(sumar_list[[i]])
# pvalues
nam <- rownames(sumar_list[[i]])
nam <- trimws(nam , which = c("both", "left", "right"))
meanSq <- sumar_list[[i]]$`Mean Sq`
names(meanSq) <- nam
Df <- sumar_list[[i]]$Df
names(Df) <- nam
pval_Candies <- pf(meanSq["Candies"]/meanSq["Candies:Judges"],
df1=Df["Candies"], df2=Df["Candies:Judges"],
lower.tail = FALSE)
pval_Judges <- pf(meanSq["Judges"]/meanSq["Residuals"],
df1=Df["Judges"], df2=Df["Residuals"],
lower.tail = FALSE)
pval_CA <- pf(meanSq["Candies:Judges"]/meanSq["Residuals"],
df1=Df["Candies:Judges"], df2=Df["Residuals"],
lower.tail = FALSE)
pval <- c(Candies=pval_Candies,
Judges=pval_Judges,
CA=pval_CA)
pvalues[i,] <- pval
}## PC 1 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 31470.6 7867.7 780.1762 < 2e-16 ***
## Judges 10 224.9 22.5 2.2304 0.02089 *
## Candies:Judges 40 707.7 17.7 1.7545 0.01158 *
## Residuals 110 1109.3 10.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC 2 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 1573.8 393.46 33.1604 < 2.2e-16 ***
## Judges 10 278.3 27.83 2.3455 0.0150274 *
## Candies:Judges 40 1053.3 26.33 2.2193 0.0005888 ***
## Residuals 110 1305.2 11.87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC 3 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 307.12 76.780 7.646 1.790e-05 ***
## Judges 10 1006.62 100.662 10.024 8.574e-12 ***
## Candies:Judges 40 484.02 12.101 1.205 0.2229
## Residuals 110 1104.61 10.042
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC 4 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 29.62 7.4041 1.0131 0.4039
## Judges 10 108.27 10.8275 1.4816 0.1557
## Candies:Judges 40 357.43 8.9357 1.2227 0.2062
## Residuals 110 803.89 7.3081
## PC 5 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 13.23 3.3069 0.7059 0.58958
## Judges 10 62.48 6.2485 1.3338 0.22151
## Candies:Judges 40 296.16 7.4040 1.5805 0.03254 *
## Residuals 110 515.32 4.6847
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC 6 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 10.99 2.7468 0.5749 0.681450
## Judges 10 133.09 13.3088 2.7852 0.004158 **
## Candies:Judges 40 176.20 4.4050 0.9219 0.605543
## Residuals 110 525.62 4.7784
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC 7 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 7.00 1.7500 0.7524 0.5585229
## Judges 10 102.49 10.2492 4.4064 3.29e-05 ***
## Candies:Judges 40 208.38 5.2094 2.2397 0.0005144 ***
## Residuals 110 255.86 2.3260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## PC 8 Df Sum Sq Mean Sq F value Pr(>F)
## Candies 4 3.594 0.8986 0.3141 0.8680
## Judges 10 32.565 3.2565 1.1383 0.3407
## Candies:Judges 40 135.848 3.3962 1.1871 0.2408
## Residuals 110 314.693 2.8608p-values
| Candies | Judges | Candies:Judges |
|---|---|---|
| 1.442e-32 | 0.02089 | 0.01158 |
| 1.495e-07 | 0.01503 | 0.0005888 |
| 0.0004731 | 8.574e-12 | 0.2229 |
| 0.5149 | 0.1557 | 0.2062 |
| 0.7742 | 0.2215 | 0.03254 |
| 0.6484 | 0.004158 | 0.6055 |
| 0.8521 | 3.29e-05 | 0.0005144 |
| 0.8989 | 0.3407 | 0.2408 |
pvalues <= 0.05
## Candies Judges Candies:Judges
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE FALSE
## [4,] FALSE FALSE FALSE
## [5,] FALSE FALSE TRUE
## [6,] FALSE TRUE FALSE
## [7,] FALSE TRUE TRUE
## [8,] FALSE FALSE FALSEBonferroni corrected p-values
pval_corrected <- t(round(pvalues*(nPC*3),4))
pval_corrected[pval_corrected>1]=1
pander(pval_corrected)| Candies | 0 | 0 | 0.0114 | 1 | 1 | 1 | 1 | 1 |
| Judges | 0.5014 | 0.3607 | 0 | 1 | 1 | 0.0998 | 8e-04 | 1 |
| Candies:Judges | 0.278 | 0.0141 | 1 | 1 | 0.781 | 1 | 0.0123 | 1 |
corrected pvalues <= 0.05
## Candies Judges Candies:Judges
## [1,] TRUE FALSE FALSE
## [2,] TRUE FALSE TRUE
## [3,] TRUE TRUE FALSE
## [4,] FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE
## [6,] FALSE FALSE FALSE
## [7,] FALSE TRUE TRUE
## [8,] FALSE FALSE FALSE6 GLM decomposition with fixed effects
6.1 Model definition
6.2 GLM decomposition
# Decomposition of Y into the effect matrices by GLM
resGLM <- matrixDecomposition(formula,outcomes,design)
modelMatrix <- sapply(resGLM$modelMatrixByEffect, function(x) x)
nparam <- sum(sapply(modelMatrix, function(x) dim(x)[2]))-1
# Building of the list of pure effects matrices
EffectMatGLM <- resGLM$effectMatrices[-1] # minus intercept
res <- vector(mode = "list")
res[[1]] <- resGLM$residuals
EffectMatGLM <- c(EffectMatGLM, residuals=res) # plus residuals
pander("names(EffectMatGLM)")names(EffectMatGLM)
## [1] "Candies" "Judges" "Candies:Judges" "residuals"6.3 scores/loadings generic graphs
listgraphs <- list()
varASCA <- list()
for(i in 1:length(ModelTerms)) {
# PCA on the pure effect matrices
ascaSVD = SVDforPCA(EffectMatGLM[[i]])
ascaSVD$scores=round(ascaSVD$scores,5)
varASCA[[i]] <- ascaSVD$var
# Scores plot
a = DrawScores(ascaSVD, type.obj = "PCA", drawNames = TRUE,
createWindow = F,
main = paste0(ModelTerms[i],"score plot - ASCA-GLM"),
color = as.factor(Judges),
pch = as.factor(Candies), axes = c(1, 2),size=2.5)
# Loadings plot
b = ScatterPlot(ascaSVD$loadings[,1],
ascaSVD$loadings[,2],
points_labs = rownames(ascaSVD$loadings),
cex.lab = 3)+
geom_vline(xintercept = 0, lwd=0.1) +
geom_hline(yintercept = 0, lwd=0.1)+
labs(title = paste0(ModelTerms[i],"- loadings plot - ASCA-GLM"),
x=paste0("PC1 (", round(ascaSVD$var[1],2), "%)"),
y=paste0("PC2 (", round(ascaSVD$var[2],2), "%)"))
listgraphs[[paste0(ModelTerms[i],"_scores")]] <- a
listgraphs[[paste0(ModelTerms[i],"_loadings")]] <- b
}
listgraphs## $Candies_scores
##
## $Candies_loadings
##
## $Judges_scores
##
## $Judges_loadings
##
## $`Candies:Judges_scores`
##
## $`Candies:Judges_loadings`
##
## $residuals_scores
##
## $residuals_loadings
7 Parallel mixed modelling
###################################################
# Parallel mixed modelling
###################################################
### Coding the interaction -------------------
CandiesJudges <- model.matrix(~ Candies:Judges-1,
data = design)
# number of parameters for the interaction
num <- 1:dim(CandiesJudges)[2]
for (i in 1:dim(CandiesJudges)[2]){
CandiesJudges[which(CandiesJudges[,i]!=0),i] = i
}
CandiesJudges <- rowSums(CandiesJudges)
CandiesJudges <- as.factor(CandiesJudges)
# New design with the coded interaction
designInter <- cbind(design,CandiesJudges)
### parallel LMM
#######################
form <- "~ Candies + (1|Judges) + (1|CandiesJudges)"
REML <- TRUE
# run parlmer -------------------
# Modification of the parlmer function (=> into parlmer_interaction)
# to obtain 2 by 2 orthogonal effect matrices.
res.parlmer <- parlmer_interaction(design=designInter, outcomes, form, REML)
MM_full <- res.parlmer
pander("summary for the first response")summary for the first response
## Linear mixed model fit by REML ['lmerMod']
##
## REML criterion at convergence: 876.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7692 -0.5925 -0.0261 0.4895 4.5514
##
## Random effects:
## Groups Name Variance Std.Dev.
## CandiesJudges (Intercept) 1.1056 1.0515
## Judges (Intercept) 0.7998 0.8943
## Residual 10.4944 3.2395
## Number of obs: 165, groups: CandiesJudges, 55; Judges, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.154e-15 3.692e-01 0.00
## Candies1 -1.684e+01 5.958e-01 -28.27
## Candies2 1.167e+01 5.958e-01 19.59
## Candies3 1.036e+01 5.958e-01 17.40
## Candies4 1.177e+01 5.958e-01 19.77
##
## Correlation of Fixed Effects:
## (Intr) Cands1 Cands2 Cands3
## Candies1 0.000
## Candies2 0.000 -0.179
## Candies3 0.000 -0.179 -0.179
## Candies4 0.000 -0.179 -0.179 -0.179# save the results -------------------
RanModMatlist <- res.parlmer$RanModMatlist
# head(RanModMatlist$`1 | CandiesJudges`)
# head(RanModMatlist$`1 | Judges`)
FixedModMatlist <- res.parlmer$FixedModMatlist
# head(FixedModMatlist$Candies)
# head(FixedModMatlist$`(Intercept)`)
# Residuals sd error
Res_std_error_PC <- sapply(MM_full$merMod_obj, sigma)
# Residuals
Residuals_PC <- sapply(MM_full$merMod_obj, residuals)
### ranef_PC: Extract the modes of the random effects
ranef_PC <- sapply(MM_full$merMod_obj, function(x) unlist(ranef(x, condVar=FALSE)))
# recover accurate rownames and colnames of ranef_PC
list_rownam <- lapply(ranef(MM_full$merMod_obj$PC1, condVar=FALSE), rownames)
colnam <- paste0(names(ranef(MM_full$merMod_obj$PC1, condVar=FALSE)),
lapply(ranef(MM_full$merMod_obj$PC1, condVar=FALSE), rownames))
colnam <- rownames(ranef_PC) #(?)
rownames(ranef_PC) <- colnam
rownames(ranef_PC) <- gsub("..(Intercept))", "", rownames(ranef_PC))
### fixef: Extract fixed-effects estimates
fixef_PC <- sapply(MM_full$merMod_obj, fixef)
### all fixed estimates and random predictions
cof_PC <- rbind(fixef_PC, ranef_PC)
### Extract Variance and Correlation Components
varcor_random_full <- sapply(MM_full$merMod_obj, VarCorr) # Std.Dev.
# varcor_random_full[,1]
# dat <- as.numeric(rbind(varcor_random_full))
# varcor_random_full_mat <- matrix(dat, nrow=2,
# dimnames = dimnames(varcor_random_full))
#
### varcor_fixed
varcor_fixed_full <- sapply(MM_full$merMod_obj, function(x) sqrt(diag(vcov(x)))) # Std.Dev.
rownames(varcor_fixed_full) <- rownames(vcov(MM_full$merMod_obj[[1]]))
##### Names of the effects
fixNames <- MM_full$fixNames # names of fixed effects
ranNames <- MM_full$ranNames # names of random effects8 Normality of the model residuals ?
library("car")
# pdf(file.path(out_path, "SDA_residuals.pdf"), height = 8, width = 8)
par(mfrow=c(4,4), mar=c(4,4,2,2))
for (i in 1:ncol(Residuals_PC)){
qqPlot(Residuals_PC[,i], main = paste0("qqplot - PC",i),
ylab="Residuals quantiles")
d = density(Residuals_PC[,i])
hist(Residuals_PC[,i], main = paste0("Histogram - PC",i), xlab="")
}
# dev.off()
9 Effect matrix computation
###################################################
# Effect matrix
###################################################
## Computation
# fixed effects + intercept -----------------------------
dim1FixedModMad <- sapply(FixedModMatlist, function(x) dim(x)[2])
names_FixedEffects <- names(FixedModMatlist)
shortFixedNames <- gsub("[^A-z]", "", names_FixedEffects)
Xmat <- do.call(cbind, FixedModMatlist)
# colnames(Xmat) %in% rownames(fixef_PC)
Xmat <- Xmat[,rownames(fixef_PC)] # reorder colnames of Xmat
index <- cumsum(dim1FixedModMad)
k <- 1
Mfix <- vector("list", length=length(shortFixedNames))
Mfix_PC <- vector("list", length=length(shortFixedNames))
names(Mfix) <- names(Mfix_PC) <- shortFixedNames
for (i in 1:length(shortFixedNames)){
XMfix = Xmat
XMfix[,-c(k:index[i])] = 0
Mfix_PC[[i]] = XMfix %*% fixef_PC
Mfix[[i]] <- Mfix_PC[[i]]%*%t(spectra_PCA_loadings)
k <- index[i] + 1
}
M0 <- Mfix[[1]]
Mfix <- Mfix[-1]
Mfix_PC <- Mfix_PC[-1]
# random effects -----------------------------
dim1RandModMad <- sapply(RanModMatlist, function(x) dim(x)[2])
names_randomEffects <- names(RanModMatlist)
shortRandNames <- gsub("[^A-z]", "", names_randomEffects)
# rbind(rownames(RanModMatlist$`1 | Judges`),rownames(RanModMatlist$`1 | CandiesJudges`)) # ok
Zmat <- do.call(cbind, RanModMatlist)
colnames(Zmat) <- paste0(rep(shortRandNames, dim1RandModMad),colnames(Zmat))
# colnames(Zmat) %in% rownames(ranef_PC)
index <- cumsum(dim1RandModMad)
k <- 1
Mrand <- vector("list", length=length(ranNames))
Mrand_PC <- vector("list", length=length(ranNames))
names(Mrand) <- names(Mrand_PC) <- ranNames
for (i in 1:length(shortRandNames)){
XMrand = Zmat
XMrand[,-c(k:index[i])] = 0
Mrand_PC[[i]] = XMrand %*% ranef_PC
Mrand[[i]] <- Mrand_PC[[i]]%*%t(spectra_PCA_loadings)
k <- index[i] + 1
}
# Residuals -----------------------------
Mres_PC <- Residuals_PC
Mres <- Mres_PC%*%t(spectra_PCA_loadings)
Mres_full <- Mres10 Are all the estimated effect matrices orthogonal ?
## t(Mrand$CandiesJudges)%*%Mrand$Judges| Â | Transp | Acid | Sweet | Raspb | Sugar |
|---|---|---|---|---|---|
| Transp | 1.582e-15 | -1.11e-16 | -9.992e-15 | -1.232e-14 | 1.277e-15 |
| Acid | 5.794e-15 | 9.176e-14 | -1.948e-14 | 2.168e-14 | 9.284e-15 |
| Sweet | 4.108e-15 | 4.441e-16 | -2.043e-14 | -3.109e-15 | 4.885e-15 |
| Raspb | 3.886e-16 | 8.438e-15 | 1.288e-14 | 1.621e-14 | -3.886e-15 |
| Sugar | -1.429e-15 | 4.108e-15 | 9.992e-15 | -1.01e-14 | -8.826e-15 |
| Bites | 8.743e-16 | -7.022e-15 | -3.22e-15 | -4.302e-16 | 2.567e-16 |
| Hard | -3.102e-15 | 9.714e-15 | 4.774e-15 | 1.221e-15 | -6.106e-16 |
| Elastic | 1.443e-15 | -7.327e-15 | 7.438e-15 | 1.332e-15 | 4.718e-16 |
| Sticky | 1.36e-15 | -4.441e-16 | -3.109e-15 | -2.887e-15 | -3.331e-16 |
| Â | Bites | Hard | Elastic | Sticky |
|---|---|---|---|---|
| Transp | -5.551e-16 | -5.551e-15 | -1.998e-15 | -3.331e-16 |
| Acid | -4.385e-15 | -1.638e-15 | 1.513e-15 | -7.216e-16 |
| Sweet | 1.954e-14 | -4.441e-16 | -2.22e-16 | -1.288e-14 |
| Raspb | 2.22e-15 | -4.885e-15 | -1.665e-15 | 7.55e-15 |
| Sugar | 3.331e-16 | -3.22e-15 | -1.221e-15 | 2.109e-15 |
| Bites | 2.914e-15 | 3.025e-15 | -1.735e-16 | 2.359e-15 |
| Hard | -5.607e-15 | -5.551e-17 | 1.693e-15 | -2.387e-15 |
| Elastic | -4.774e-15 | 4.108e-15 | -1.055e-15 | -2.554e-15 |
| Sticky | 6.883e-15 | 6.661e-16 | 2.554e-15 | 5.107e-15 |
# unique(Ma[,1])
# unique(Mca[,1])
#
# colSums(Ma)
# colSums(Mca)
cat("t(Mfix$Candies)%*%(Mrand$CandiesJudges)")## t(Mfix$Candies)%*%(Mrand$CandiesJudges)| Â | Transp | Acid | Sweet | Raspb | Sugar |
|---|---|---|---|---|---|
| Transp | 8.846e-13 | 1.088e-13 | -1.3e-12 | -1.03e-12 | -7.994e-13 |
| Acid | -8.651e-13 | -6.262e-14 | 1.116e-12 | 8.917e-13 | 7.656e-13 |
| Sweet | -2.269e-13 | -5.052e-14 | 3.784e-13 | 3.046e-13 | 1.852e-13 |
| Raspb | -2.891e-13 | -2.909e-14 | 4.032e-13 | 3.411e-13 | 2.42e-13 |
| Sugar | -9.273e-13 | -1.292e-13 | 1.386e-12 | 1.073e-12 | 8.562e-13 |
| Bites | 8.278e-13 | 1.181e-13 | -1.201e-12 | -9.45e-13 | -7.301e-13 |
| Hard | 8.278e-13 | 9.992e-14 | -1.187e-12 | -9.486e-13 | -7.354e-13 |
| Elastic | 8.153e-13 | 1.137e-13 | -1.208e-12 | -9.486e-13 | -7.407e-13 |
| Sticky | 7.327e-13 | 1.035e-13 | -1.03e-12 | -8.562e-13 | -6.537e-13 |
| Â | Bites | Hard | Elastic | Sticky |
|---|---|---|---|---|
| Transp | 8.786e-13 | 8.766e-13 | 7.514e-13 | 7.496e-13 |
| Acid | -8.484e-13 | -8.695e-13 | -7.621e-13 | -7.532e-13 |
| Sweet | -2.42e-13 | -2.338e-13 | -1.779e-13 | -1.767e-13 |
| Raspb | -2.874e-13 | -2.9e-13 | -2.3e-13 | -2.549e-13 |
| Sugar | -9.213e-13 | -9.335e-13 | -7.834e-13 | -7.816e-13 |
| Bites | 8.298e-13 | 8.251e-13 | 6.883e-13 | 6.892e-13 |
| Hard | 8.096e-13 | 8.251e-13 | 6.99e-13 | 7.141e-13 |
| Elastic | 8.107e-13 | 8.162e-13 | 6.883e-13 | 7.141e-13 |
| Sticky | 7.427e-13 | 7.443e-13 | 6.102e-13 | 6.217e-13 |
11 PCA on the effect matrices
11.1 PCA on the pure effect matrices
###################################################
## PCA on the effect matrices
###################################################
## Fixed effects -------------------------------
loadingsFix <- vector(mode = "list", length = length(Mfix))
for (i in 1:length(Mfix)){
print(names(Mfix[i]))
dimnames(Mfix[[i]]) <- dimnames(raw_outcomes)
res_pca <- SVDforPCA(Mfix[[i]])
res_pca$scores <- round(res_pca$scores, 6)
barplot(res_pca$var[1:10], main = "scree plot")
print(DrawScores(res_pca, size=3, color = design[,names(Mfix[i])],
main = paste("Scores plot"),
axes=c(1,2), drawNames=FALSE))
print(DrawScores(res_pca, size=3, color = design[,names(Mfix[i])],
main = paste("Scores plot"), axes=c(3,4),
drawNames=FALSE))
load <- DrawLoadings(res_pca, type.obj = "PCA", createWindow = F,
main = paste("Loadings plot"),
axes = c(1:4), loadingstype = "s",
num.stacked = 4, nxaxis = 9,
ang = "0",
xaxis_type = "character")
load2 <- ScatterPlot(res_pca$loadings[,1],
res_pca$loadings[,2],
points_labs = rownames(res_pca$loadings),
cex.lab = 5)
load3 <- ScatterPlot(res_pca$loadings[,3],
res_pca$loadings[,4],
points_labs = rownames(res_pca$loadings),
cex.lab = 5)
print(load[[1]])
print(load2)
print(load3)
loadingsFix[[i]] <- load
}## [1] "Candies"
## Random effects -------------------------------
loadingsRand <- vector(mode = "list", length = length(Mrand))
for (i in 1:length(Mrand)){
print(names(Mrand[i]))
dimnames(Mrand[[i]]) <- dimnames(raw_outcomes)
res_pca <- SVDforPCA(Mrand[[i]])
barplot(res_pca$var[1:10], main = "scree plot")
res_pca$scores <- round(res_pca$scores, 6)
print(DrawScores(res_pca, size=3, color = design$Judges,
pch=design$Candies,
main = paste("Scores plot"),
axes=c(1,2), drawNames=FALSE))
print(DrawScores(res_pca, size=3, color = design$Judges,
pch=design$Candies,
main = paste("Scores plot"), axes=c(3,4),
drawNames=FALSE))
load <- DrawLoadings(res_pca, axes=c(1:4), type.obj = "PCA",
xaxis_type = "character",
loadingstype = "s", xlab = "ppm", nxaxis = 9,
main = paste("Loadings plot:"))
load2 <- ScatterPlot(res_pca$loadings[,1],
res_pca$loadings[,2],
points_labs = rownames(res_pca$loadings),
cex.lab = 5)
load3 <- ScatterPlot(res_pca$loadings[,3],
res_pca$loadings[,4],
points_labs = rownames(res_pca$loadings),
cex.lab = 5)
print(load[[1]])
print(load2)
print(load3)
loadingsRand[[i]] <- load
}## [1] "CandiesJudges"
## [1] "Judges"
11.1.1 Scree plots and pure loadings
Note that only the loadings are backtransformed
# ========================
# Pure loadings ==========
# ========================
# Candies
#####################
res_pca <- SVDforPCA(Mfix_PC$Candies)
colnames(Mfix_PC$Candies) <- paste0("Var",colnames(Mfix_PC$Candies))
df <- data.frame(PC = as.character(1:8), var = res_pca$var)
screeplotCandies <- ggplot(df, aes(y=0,yend=var,x = PC,
xend=PC))+ geom_segment()+
labs(title= "Scree plot",
x = "PC", y="% var") + theme_classic()
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadCandies <- ScatterPlot(backtransf_load[,1],
backtransf_load[,2],
points_labs = rownames(backtransf_load),
cex.lab = 3)+
geom_vline(xintercept = 0, lwd=0.1) +
geom_hline(yintercept = 0, lwd=0.1)+
coord_cartesian(xlim = c(-0.5, 0.5), ylim = c(-1, 1))+
labs(title = "PCA loadings", x=paste0("PC1 (", round(res_pca$var[1],2), "%)"),
y=paste0("PC2 (", round(res_pca$var[2],2), "%)"))
## Random effects -------------------------------
# Judges ++++++++
res_pca <- SVDforPCA(Mrand_PC$Judges)
df <- data.frame(PC = as.character(1:length(res_pca$var)),
var = res_pca$var)
screeplotJudges <- ggplot(df, aes(y=0,yend=var,x=PC,
xend=PC))+ geom_segment()+
labs(title= "Scree plot",
x = "PC", y="% var")+theme_classic()
# screeplot
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadJudges <- ScatterPlot(backtransf_load[,1],
backtransf_load[,2],
points_labs = rownames(backtransf_load),
cex.lab = 3)+
geom_vline(xintercept = 0, lwd=0.1) +
geom_hline(yintercept = 0, lwd=0.1)+
coord_cartesian(xlim = c(-0.8, 0.8), ylim = c(-0.6, 0.6))+
labs(title = "PCA loadings", x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)"))
# CandiesJudges ++++++++
res_pca <- SVDforPCA(Mrand_PC$CandiesJudges)
df <- data.frame(PC = as.character(1:length(res_pca$var)),
var = res_pca$var)
screeplotCandiesJudges <- ggplot(df, aes(y=0,yend=var,x=PC,
xend=PC))+ geom_segment()+
labs(title= "Scree plot",
x = "PC", y="% var")+theme_classic()
# screeplot
backtransf_load <- -1*t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadCandiesJudges <- ScatterPlot(backtransf_load[,1],
backtransf_load[,2],
points_labs = rownames(backtransf_load),
cex.lab = 3) +
geom_vline(xintercept = 0, lwd=0.1) +
geom_hline(yintercept = 0, lwd=0.1)+
coord_cartesian(xlim = c(-0.8, 0.8), ylim = c(-0.6, 0.6)) +
labs(title = "PCA loadings", x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)"))11.2 PCA on the augmented effect matrices (ASCA-E)
11.2.1 addSegments fun
addSegments <- function(group, data, pch=16, main = NULL,
col = rainbow(n = length(unique(group))),
...) {
group <- as.factor(group)
# col <- rainbow(n = length(unique(group)))
plot(data$x, data$y, col=col[group], pch=pch, main=main, ...)
group <- as.factor(group)
xcent <- tapply(data[,1], group, FUN=mean)
ycent <- tapply(data[,2], group, FUN=mean)
centers <- data.frame(xcent=xcent, ycent=ycent)
mapply(FUN = points, x = centers$xcent, y = centers$ycent,
col = col, MoreArgs = list(pch=20, cex=0.7))
submatrices <- split(x=data, f=group)
for (i in 1:nlevels(group)){
mapply(FUN = segments, x1 = submatrices[[i]]$x,
y1 = submatrices[[i]]$y, col = col[i],
MoreArgs = list(x0 = centers$xcent[i],
y0 = centers$ycent[i]))
}
}11.2.2 Augmented (ANOVA) scores plots
# ========================
# augmented scores =======
# ========================
# * Candies: Mc + Mca
# * Assessor: Ma + E
# * Candies:Assessor : Mca + E
a <- nlevels(design$Candies)
b <- nlevels(design$Judges)
nn <- table(design$Judges,design$Candies)[1,1]
# Candies
##################
# pander("Candies")
ascaSVD = SVDforPCA(Mfix_PC$Candies)
Fstat <- qf(.95, df1=(a-1), df2=((a-1)*(b-1)))
coef <- sqrt(Fstat/(b-1))
ascaSVD$scores[,1:2]=(Mfix_PC$Candies +
((Mrand_PC$CandiesJudges)*coef)) %*%
ascaSVD$loadings[,1:2]
Candies_scores <- DrawScores(ascaSVD,
main = "PCA Scores",
color = Candies, pch = Candies,
drawNames = FALSE, drawPolygon = FALSE, drawEllipses = TRUE,
noLegend = FALSE, size=2)
# Judges
##################
# pander("Judges")
ascaSVD = SVDforPCA(Mrand_PC$Judges)
df1 <- (b-1)
df2 <- (a*b*(nn-1))
Fstat <- qf(.95, df1=df1, df2=df2)
coef <- sqrt(Fstat*df1/df2)
ascaSVD$scores[,1:2]=(Mrand_PC$Judges+(Mres_PC*coef))%*%
ascaSVD$loadings[,1:2]
Judges_scores <- DrawScores(ascaSVD, type.obj = "PCA",
drawNames = F, createWindow = F,
main = "PCA Scores",
color = Judges,
pch = Judges, axes = c(1, 2),size=2, drawPolygon = FALSE,
drawEllipses = TRUE)+
theme(legend.text=element_text(size=10),
legend.key.height=unit(0.7,"line"))
# CandiesJudges
##################
# pander("CandiesJudges")
ascaSVD = SVDforPCA(Mrand_PC$CandiesJudges)
df1 <- (a-1)*(b-1)
df2 <- (a*b*(nn-1))
Fstat <- qf(.95, df1=df1, df2=df2)
coef <- sqrt(Fstat*df1/df2)
ascaSVD$scores[,1:2]=-1*(Mrand_PC$CandiesJudges+(Mres_PC*coef))%*%
ascaSVD$loadings[,1:2]
CJ <- as.factor(paste0(design$Candies, design$Judges))
index <- c(145,147,27, 109, 110)
data <- ascaSVD$scores[,1:2]
group <- as.factor(CJ)
xcent <- tapply(data[,1], group, FUN=mean)
ycent <- tapply(data[,2], group, FUN=mean)
xc <- yc <- c()
for (i in 1:length(group)){
xc[i] <- xcent[group[i]]
yc[i] <- ycent[group[i]]
}
df <- data.frame(data, xcent=xc, ycent = yc, CJ=CJ,
Judges =design$Judges, Candies=design$Candies)
# Eigenvalues
eig <- ascaSVD$eigval
# Variances in percentage
variance <- eig * 100/sum(eig)
Xax=1
Yax=2
XaxName <- paste0("PC", Xax, " (", round(variance[Xax], 2),"%)")
YaxName <- paste0("PC", Yax, " (", round(variance[Yax], 2), "%)")
xlab <- XaxName
ylab <- YaxName
Xlim <- c(min(ascaSVD$scores[, Xax]) * 1.4, max(ascaSVD$scores[, Xax]) * 1.4)
Ylim <- c(min(ascaSVD$scores[, Yax]) * 1.4, max(ascaSVD$scores[, Yax]) * 1.4)
main = "PCA Scores"
plots <- ggplot(df,aes(x=xc,y=yc))+
geom_point(df, mapping=aes(x=PC1,y=PC2, shape=Candies, color=Judges),size=2) +
geom_segment(aes(yend=PC2,xend=PC1,color=Judges,group=CJ)) +
ggplot2::xlim(Xlim) + ggplot2::ylim(Ylim) +
scale_shape_manual(values=seq(0,26), name = "Candies")
plots <- plots + ggplot2::labs(title = main, x = xlab, y = ylab) +
ggplot2::geom_vline(xintercept = 0,size = 0.1) +
ggplot2::geom_hline(yintercept = 0, size = 0.1) + ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major =
ggplot2::element_line(color = "gray60", size = 0.2),
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "gray98"))
plots <- plots + theme(legend.text=element_text(size=10),
legend.key.height=unit(0.6,"line")) +
annotate("text", y = (ascaSVD$scores[index,2] +
1.7*c(-1 ,1,-1,-1, -1)),
x = ascaSVD$scores[index,1],
label = row.names(outcomes[index,1:2]))
CA_scores <- plots
mat <- data.frame(x=ascaSVD$scores[,1],y=ascaSVD$scores[,2])
group <- as.factor(paste0(Candies,Judges))
# table(group)
# pdf(file.path(fig_path,"SDA_ScoresCJEM_linkbyCJ.pdf"))
# addSegments(group=group, data=mat, main = "Scores plot C*J effect matrix \n color by CandiesJudges", ylab = "PC2", xlab = "PC1")
# abline(h=0, v=0, lty=2, lwd =1,col="black")
# dev.off()
# Residuals ================
res_pca <- SVDforPCA(Mres_PC)
# colnames(Mfix_PC$Candies) <- paste0("Var",colnames(Mfix_PC$Candies))
df <- data.frame(PC = as.character(1:8), var = res_pca$var)
# df <- df[1:7,]
screeplot_resid <- ggplot(df, aes(y=0,yend=var,x = PC,
xend=PC))+ geom_segment()+
labs(title= "Scree plot",
x = "PC", y="% var") + theme_classic()
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadings_resid <- ScatterPlot(backtransf_load[,1],
backtransf_load[,2],
points_labs = rownames(backtransf_load),
cex.lab = 3)+
geom_vline(xintercept = 0, lwd=0.1) +
geom_hline(yintercept = 0, lwd=0.1)+
labs(title = "PCA loadings",
x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)")) +
coord_cartesian(xlim = c(-0.85, 0.85), ylim = c(-0.5, 0.5))
index <- c(10, 27)
df <- data.frame(PC = as.character(1:length(res_pca$var)),
var = res_pca$var)
df <- df[1:7,]
scores_resid <- DrawScores(res_pca, drawNames = FALSE,
color = Candies,
pch = Judges, main ="PCA scores",
size = 2) +
coord_cartesian(xlim = c(-16, 16), ylim = c(-18, 18))+
theme(legend.text=element_text(size=10), legend.key.height=unit(0.7,"line"))+
annotate("text", y = (res_pca$scores[index,2] +
1.7*c(-1 , -1)),
x = res_pca$scores[index,1],
label = rownames(res_pca$scores[index,1:2]))11.2.3 Augmented (ED) scores plots
Effective dimensions
EDc <- rep((nlevels(design$Candies)-1),ncol(outcomes))
EDj <- ED_sensory["ED_judges",]
EDcj <- ED_sensory["ED_cj",]
EDe <- n - colSums(ED_sensory)
# ========================
# augmented scores =======
# ========================
# * Candies: Mc + Mca
# * Assessor: Ma + E
# * Candies:Assessor : Mca + E
# Candies
##################
# pander("Candies")
ascaSVD = SVDforPCA(Mfix_PC$Candies)
df1 <- EDc
df2= EDcj
Fstat <- qf(.95, df1=df1 ,df2= df2)
coef <- sqrt((Fstat*df1)/df2)
mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
mat[,i] <- Mrand_PC$CandiesJudges[,i]*coef[i]
}
ascaSVD$scores[,1:2]= (Mfix_PC$Candies + mat) %*%
ascaSVD$loadings[,1:2]
Candies_scores <- DrawScores(ascaSVD,
main = "PCA Scores",
color = Candies, pch = Candies,
drawNames = FALSE, drawPolygon = FALSE, drawEllipses = TRUE,
noLegend = FALSE, size=2)
# Judges
##################
# pander("Judges")
ascaSVD = SVDforPCA(Mrand_PC$Judges)
df1 <- EDj
df2 <- EDe
Fstat <- qf(.95, df1=df1, df2=df2)
coef <- sqrt(Fstat*df1/df2)
mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
mat[,i] <- Mres_PC[,i]*coef[i]
}
ascaSVD$scores[,1:2]=(Mrand_PC$Judges + mat)%*%
ascaSVD$loadings[,1:2]
Judges_scores <- DrawScores(ascaSVD, type.obj = "PCA",
drawNames = F, createWindow = F,
main = "PCA Scores",
color = Judges,
pch = Judges, axes = c(1, 2),size=2, drawPolygon = FALSE,
drawEllipses = TRUE)+
theme(legend.text=element_text(size=10),
legend.key.height=unit(0.7,"line"))
# CandiesJudges
##################
# pander("CandiesJudges")
ascaSVD = SVDforPCA(Mrand_PC$CandiesJudges)
df1 <- EDcj
df2 <- EDe
Fstat <- qf(.95, df1=df1, df2=df2)
coef <- sqrt(Fstat*df1/df2)
mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
mat[,i] <- Mres_PC[,i]*coef[i]
}
ascaSVD$scores[,1:2]=-1*(Mrand_PC$CandiesJudges + mat)%*%
ascaSVD$loadings[,1:2]
CJ <- as.factor(paste0(design$Candies, design$Judges))
index <- c(145,147,27, 109, 110)
data <- ascaSVD$scores[,1:2]
group <- as.factor(CJ)
xcent <- tapply(data[,1], group, FUN=mean)
ycent <- tapply(data[,2], group, FUN=mean)
xc <- yc <- c()
for (i in 1:length(group)){
xc[i] <- xcent[group[i]]
yc[i] <- ycent[group[i]]
}
df <- data.frame(data, xcent=xc, ycent = yc, CJ=CJ,
Judges =design$Judges, Candies=design$Candies)
# Eigenvalues
eig <- ascaSVD$eigval
# Variances in percentage
variance <- eig * 100/sum(eig)
Xax=1
Yax=2
XaxName <- paste0("PC", Xax, " (", round(variance[Xax], 2),"%)")
YaxName <- paste0("PC", Yax, " (", round(variance[Yax], 2), "%)")
xlab <- XaxName
ylab <- YaxName
Xlim <- c(min(ascaSVD$scores[, Xax]) * 1.4, max(ascaSVD$scores[, Xax]) * 1.4)
Ylim <- c(min(ascaSVD$scores[, Yax]) * 1.4, max(ascaSVD$scores[, Yax]) * 1.4)
main = "PCA Scores"
plots <- ggplot(df,aes(x=xc,y=yc))+
geom_point(df, mapping=aes(x=PC1,y=PC2, shape=Candies, color=Judges),size=2) +
geom_segment(aes(yend=PC2,xend=PC1,color=Judges,group=CJ)) +
ggplot2::xlim(Xlim) + ggplot2::ylim(Ylim) +
scale_shape_manual(values=seq(0,26), name = "Candies")
plots <- plots + ggplot2::labs(title = main, x = xlab, y = ylab) +
ggplot2::geom_vline(xintercept = 0,size = 0.1) +
ggplot2::geom_hline(yintercept = 0, size = 0.1) + ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major =
ggplot2::element_line(color = "gray60", size = 0.2),
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_rect(fill = "gray98"))
plots <- plots + theme(legend.text=element_text(size=10),
legend.key.height=unit(0.6,"line"))
CA_scores <- plots
# CA_scores
mat <- data.frame(x=ascaSVD$scores[,1],y=ascaSVD$scores[,2])
group <- as.factor(paste0(Candies,Judges))
# table(group)
# pdf(file.path(fig_path,"SDA_ScoresCJEM_linkbyCJ.pdf"))
# addSegments(group=group, data=mat, main = "Scores plot C*J effect matrix \n color by CandiesJudges", ylab = "PC2", xlab = "PC1")
# abline(h=0, v=0, lty=2, lwd =1,col="black")
# dev.off()
# ========================
# Residuals ==============
# ========================
res_pca <- SVDforPCA(Mres_PC)
# colnames(Mfix_PC$Candies) <- paste0("Var",colnames(Mfix_PC$Candies))
df <- data.frame(PC = as.character(1:8), var = res_pca$var)
# df <- df[1:7,]
screeplot_resid <- ggplot(df, aes(y=0,yend=var,x = PC,
xend=PC))+ geom_segment()+
labs(title= "Scree plot",
x = "PC", y="% var") + theme_classic()
backtransf_load <- t(t(res_pca$loadings) %*% t(spectra_PCA_loadings))
loadings_resid <- ScatterPlot(backtransf_load[,1],
backtransf_load[,2],
points_labs = rownames(backtransf_load),
cex.lab = 3)+
geom_vline(xintercept = 0, lwd=0.1) +
geom_hline(yintercept = 0, lwd=0.1)+
labs(title = "PCA loadings",
x=paste0("PC1 (", round(res_pca$var[1],2), "%)"), y=paste0("PC2 (", round(res_pca$var[2],2), "%)")) +
coord_cartesian(xlim = c(-0.85, 0.85), ylim = c(-0.5, 0.5))
index <- c(10, 27)
df <- data.frame(PC = as.character(1:length(res_pca$var)),
var = res_pca$var)
df <- df[1:7,]
scores_resid <- DrawScores(res_pca, drawNames = FALSE,
color = Candies,
pch = Judges, main ="PCA scores",
size = 2) +
coord_cartesian(xlim = c(-19, 19), ylim = c(-18, 18))+
theme(legend.text=element_text(size=10), legend.key.height=unit(0.7,"line"))+
annotate("text", y = (res_pca$scores[index,2] +
1.7*c(-1 , -1)),
x = res_pca$scores[index,1],
label = rownames(res_pca$scores[index,1:2]))11.2.4 Graphs
Scores_Loadings_EffectMat <- ggarrange(screeplotCandies,
Candies_scores,
loadCandies,
screeplotJudges,
Judges_scores,
loadJudges,
screeplotCandiesJudges,
CA_scores,
loadCandiesJudges,
labels = c("A", "B", "C", "D", "E", "F", "G", "H", "I"),
ncol = 3, nrow = 3, common.legend=TRUE, widths = c(0.8,2,2))
# Scores_Loadings_EffectMat
# ggexport(Scores_Loadings_EffectMat, filename = file.path(fig_path,"SDA_Scores_Loadings_EffectMat.pdf"),
# height = 12, width = 12)
a <- grid.arrange(screeplotCandies, Candies_scores,
loadCandies, nrow=1,widths=c(0.3, 1, 0.85),
top=textGrob("Candy effect matrix",
gp=gpar(fontsize=20,font=2)))
b <- grid.arrange(screeplotJudges, Judges_scores,
loadJudges, nrow=1,widths=c(0.3, 1, 0.85),
top=textGrob("Judge effect matrix",
gp=gpar(fontsize=20,font=2)))
c <- grid.arrange(screeplotCandiesJudges,
CA_scores, loadCandiesJudges,
nrow=1,widths = c(0.3, 1, 0.85),
top=textGrob("C*J effect matrix",
gp=gpar(fontsize=20,font=2)))
d <- grid.arrange(screeplot_resid,
scores_resid, loadings_resid,
nrow=1,widths = c(0.3, 1, 0.85),
top=textGrob("Residual effect matrix",
gp=gpar(fontsize=20,font=2)))
# Thesis chapter output ----------
p <- ggarrange(a,b,c, nrow=3, labels = c("A", "B", "C"))
# ggsave(file.path(fig_path,"SDA_Scores_Loadings_EffectMat.pdf"),
# plot = p, height = 15, width = 14,scale = 0.65)
# journal article output ----------
plots <- ggarrange(a,b,c,d, nrow=4, labels = c("A", "B", "C", "D"))
# ggsave(file.path(fig_path,"SDA_Scores_Loadings_EffectMat.pdf"),
# plot = plots, height = 20, width = 14, scale = 0.7)11.2.5 Spotted outliers
# outliers detected
index <- which(design$Judges=="07" & design$Candies=="1")
tab_a <- raw_outcomes[index,]
rownames(tab_a) <- paste0("07_1_",1:3)
index <- which(design$Judges=="02" & design$Candies=="1")
tab_b <- raw_outcomes[index,]
rownames(tab_b) <- paste0("02_1_",1:3)
# ggexport(ggtexttable(tab_a, theme = ttheme("classic")),
# filename =
# file.path(fig_path,"SDA_Residuals_outcomes_07_1.pdf"),
# height = 2, width = 6)
# ggexport(ggtexttable(tab_b, theme = ttheme("classic")),
# filename =
# file.path(fig_path,"SDA_Residuals_outcomes_02_1.pdf"),
# height = 2, width = 6)11.2.6 Interaction graphs
11.2.6.1 Pure scores
# 1. Pure scores
#############################
res_pca <- SVDforPCA(Mrand_PC$CandiesJudges)
id <- which(!duplicated(res_pca$scores[,1]))
df <- as.data.frame(cbind(Candies=Candies[id],
Judges = Judges[id],
scores = res_pca$scores[id,]))
df$Candies <- as.factor(df$Candies)
df$Judges <- as.factor(df$Judges)
a <- ggplot(data=df, aes(x=Candies, y=PC1)) +
geom_point(aes(colour = Judges, shape = Judges), size=3) +
scale_shape_manual(values = c(0:10)) +
geom_hline(yintercept=0, linetype=3, color = "black")
b <- ggplot(data=df, aes(x=Judges, y=PC1)) +
geom_point(aes(colour = Candies, shape =Candies), size=3)+
geom_hline(yintercept=0, linetype=3, color = "black")
c <- ggplot(data=df, aes(x=Candies, y=PC2)) +
geom_point(aes(colour = Judges, shape = Judges), size=3) +
scale_shape_manual(values = c(0:10))+
geom_hline(yintercept=0, linetype=3, color = "black")
d <-ggplot(data=df, aes(x=Judges, y=PC2))+
geom_point(aes(colour = Candies, shape =Candies), size=3)+
geom_hline(yintercept=0, linetype=3, color = "black")# ggexport(ggarrange(a, c, common.legend=TRUE, legend="right"),
# filename = file.path(fig_path,"SDA_Scores_CAinteraction.pdf"),
# height = 4, width = 12)
# pdf(file = file.path(fig_path,"SDA_Scores_CAinteraction.pdf"),
# height = 4, width = 12)
grid.arrange(a, c, nrow=1,widths=c(1, 0.85),
top=textGrob("PCA scores on augmented C*J effect matrix",gp=gpar(fontsize=20,font=2)))
11.2.6.2 Augmented scores
# 2. Augmented scores
#############################
# a <- nlevels(design$Candies)
# b <- nlevels(design$Judges)
# nn <- table(design$Judges,design$Candies)[1,1]
EDc <- rep((nlevels(design$Candies)-1),ncol(outcomes))
EDj <- ED_sensory["ED_judges",]
EDcj <- ED_sensory["ED_cj",]
EDe <- n - colSums(ED_sensory)
pander("CandiesJudges")CandiesJudges
res_pca = SVDforPCA(Mrand_PC$CandiesJudges)
# df1 <- (a-1)*(b-1)
# df2 <- (a*b*(nn-1))
# Fstat <- qf(.95, df1=df1, df2=df2)
# coef <- sqrt(Fstat*df1/df2)
#
# res_pca$scores[,1:2]=-1*(Mrand_PC$CandiesJudges+(Mres_PC*coef))%*%
# res_pca$loadings[,1:2]
df1 <- EDcj
df2 <- EDe
Fstat <- qf(.95, df1=df1, df2=df2)
coef <- sqrt(Fstat*df1/df2)
mat <- matrix(NA, ncol=ncol(outcomes), nrow=nrow(outcomes))
for (i in 1:ncol(outcomes)){
mat[,i] <- Mres_PC[,i]*coef[i]
}
res_pca$scores[,1:2]=-1*(Mrand_PC$CandiesJudges + mat)%*%
res_pca$loadings[,1:2]
df <- as.data.frame(cbind(Candies=Candies,
Judges = Judges,
scores = res_pca$scores))
df$Candies <- as.factor(df$Candies)
df$Judges <- as.factor(df$Judges)
a <- ggplot(data=df, aes(x=Candies, y=PC1)) + geom_point(aes(colour = Judges, shape = Judges), size=3) +
scale_shape_manual(values = c(0:10)) +
geom_hline(yintercept=0, linetype=3, color = "black")+
theme(text=element_text(size=15),
legend.key.height = unit("0.7", units = "line"))
b <- ggplot(data=df, aes(x=Judges, y=PC1)) + geom_point(aes(colour = Candies, shape =Candies), size=3)+
geom_hline(yintercept=0, linetype=3, color = "black")+
theme(text=element_text(size=15),
legend.key.height = unit("0.7", units = "line"))
# interPC1 <- grid.arrange(a, b)
c <- ggplot(data=df, aes(x=Candies, y=PC2)) + geom_point(aes(colour = Judges, shape = Judges), size=3) +
scale_shape_manual(values = c(0:10))+
geom_hline(yintercept=0, linetype=3, color = "black")+
theme(text=element_text(size=15),
legend.key.height = unit("0.7", units = "line"))
d <-ggplot(data=df, aes(x=Judges, y=PC2))+geom_point(aes(colour = Candies, shape =Candies), size=3)+
geom_hline(yintercept=0, linetype=3, color = "black")+
theme(text=element_text(size=15),
legend.key.height = unit("0.7", units = "line"))
# interPC2 <- grid.arrange(c, d)
a
# ggexport(ggarrange(a, c, common.legend=TRUE, legend="right"),
# filename = file.path(fig_path,"SDA_Scores_CAinteraction_aug.pdf"),
# height = 4, width = 12)
p <- grid.arrange(a, c, nrow=1,widths=c(1, 1),
top=textGrob("PCA scores on augmented C*J effect matrix",
gp=gpar(fontsize=20,font=2)))
12 Percentage of explained variance
12.1 Method from Nakagawa and Schielzeth (2012)
####################################################
# Percentage of explained variance
####################################################
## Method from Nakagawa and Schielzeth (2012)
# Random effects -----------------------------------
sigma2_res = Res_std_error_PC^2 # Residual
varcor_random_full <- as.data.frame(varcor_random_full)
Var_Mrand <- rbind(varcor_random_full, sigma2_res=sigma2_res) # only random effects
Var_Mrand <- data.matrix(Var_Mrand)
# fixed effect -----------------------------------
Var_Mfix <- c()
for (i in 1:length(fixNames)){
# variance of parameters values (population)
Var_Mfix <- rbind(Var_Mfix, (apply(Mfix_PC[[i]], 2, var) *(n - 1) / n))
}
# all together -----------------------------------
rownames(Var_Mfix) <- fixNames
rownames(Var_Mrand) <- c(ranNames, "Residuals")
var_comp <- rbind(Var_Mfix, Var_Mrand)
var_comp <- var_comp[c(1,3,2,4),]
# log of var comp +++++
log_var_comp <- t(log1p(var_comp)) # log(x+1)
log_var_comp <- cbind(id=rownames(log_var_comp), log_var_comp)
log_var_comp <- as.data.frame(log_var_comp)
log_var_comp <- melt(log_var_comp, id=c("id"))
log_var_comp$value <- as.numeric(log_var_comp$value)
names(log_var_comp) <- c("PC", "Effect", "Variance")
sum_var_comp <- rowSums(var_comp)
# Percent var expl by each effect -----
var_components_abs <- sum_var_comp
var_components <- var_components_abs*100/sum(var_components_abs)
names(var_components) <- rownames(var_comp)
# pdf(file.path(fig.path,"variance_components.pdf"), height = 3, width = 3, pointsize = 12)
par(mar=c(0.5,2,3,0.5))
barplot(var_components, main="Variance components \n percentage",xaxt="n",las=2, col=c(darkblue,turquoise,violetred , limegreen, gray67), border = NA,
legend = names(var_components), args.legend = list(x="topright",
inset=c(0,0),box.lty=0,cex = 1, y.intersp = 0.8))
# dev.off()
table_var <- rbind(var_components_abs, var_components)
rownames(table_var) <- c("Sum variance for all responses", "percentage of variation")
colnames(table_var) <- rownames(var_comp)
pander(table_var)| Â | Candies | Judges | CandiesJudges | Residuals |
|---|---|---|---|---|
| Sum variance for all responses | 202.5 | 9.442 | 5.644 | 53.24 |
| percentage of variation | 74.77 | 3.486 | 2.084 | 19.66 |
12.2 Graph
log_var_comp$PC <- sub("Var", "", log_var_comp$PC)
p <- ggplot(log_var_comp, aes(Effect, Variance, group = PC))+
ggtitle("SDA - Log of variance components")
p <- p + geom_point(aes(colour = PC))+
geom_line(aes(colour = PC, linetype=PC),size=0.5)+
theme_classic()+
theme(legend.key.width = unit(0.8,"cm"),
plot.title = element_text(size = 17),
axis.title.y = element_text(size = 17),
axis.title.x = element_text(size = 17),
axis.text = element_text(size = 13),
legend.text = element_text(size = 12),
legend.title = element_text(size = 17)) +
ylab(label = "log(variance)")
# Thesis chapter output
# ggexport(p, filename = file.path(fig_path,"SDA_variance_components.pdf"),
# height = 5, width = 6.5, pointsize=20)
tab <- data.frame(Effect= names(var_components),
pcvar = round(var_components,2))
var_comp.table <- ggtexttable(tab, cols = c("Effect", "Global var (%)"),
rows = NULL,
theme = ttheme("classic",base_size = 10))
p <- p + annotation_custom(ggplotGrob(var_comp.table),
xmin = 6, ymin = 3,
xmax = 0) +
theme(legend.title = element_blank())
# journal article output
# ggexport(p, filename = file.path(fig_path,"SDA_variance_components.pdf"),
# height = 5, width = 6.5, pointsize=19)
p
13 Bootstrap
13.1 Set up
# set up of the bootstrap
set.seed(2018)
nsim = 2000 # number of simulations
# name of the output
name_RData <- "bootstrap_Sensory_Data.RData"
# Set up -------------------
# formulas without the effet to test
null_formulas <- list(Candies = "~ (1|Judges) + (1|CandiesJudges)",
Judges = "~ Candies + (1|CandiesJudges)",
CandiesJudges = "~Candies + (1|Judges)")
null_effects <- names(null_formulas)
REML <- c(FALSE, TRUE, TRUE)
names(REML) <- null_effects13.2 log-likelihood Ratio computation
13.2.1 True log-likelihood Ratio statistics
##################################################
# True log-likelihood Ratio statistics #
##################################################
# full model: MM_full
######################
# REML +++++
loglik_PC_full_REML <- sapply(MM_full$merMod_obj, logLik, REML=T)
# ML +++++
loglik_PC_full_ML <- sapply(MM_full$merMod_obj, logLik, REML=F)
loglik_PC_full <- matrix(NA, ncol = nPC,
nrow=length(null_formulas), byrow = TRUE)
for (i in 1:length(REML)){
if (REML[i]==TRUE){
loglik_PC_full[i,] <- loglik_PC_full_REML
}else {loglik_PC_full[i,] <- loglik_PC_full_ML}
}
### Restricted models
######################
res.parlmer_NULL <- vector("list", length = length(null_formulas))
names(res.parlmer_NULL) <- names(REML) <- names(null_formulas)
for (i in 1:length(null_formulas)) {
# run parlmer
res.parlmer_NULL[[i]] <- parlmer(designInter,
outcomes, null_formulas[[i]], REML=REML[i])
}
# Save the results ------------
MM_PC_null <- lapply(res.parlmer_NULL, function(x) x[["merMod_obj"]])
#### Randvarnames, Fixvarnames
ranNames <- sapply(res.parlmer_NULL, function(x) x[["ranNames"]])
Fixvarnames <- sapply(res.parlmer_NULL, function(x) x[["fixNames"]])
varcor_random <- vector(mode = "list", length = length(null_formulas))
fixef_PC <- vector(mode = "list", length = length(null_formulas))
modmat_fixed <- vector(mode = "list", length = length(null_formulas))
M0 <- vector(mode = "list", length = length(null_formulas))
for (i in 1:length(null_formulas)){
varcor_random[[i]] <- sapply(MM_PC_null[[i]], function(x)
as.data.frame(VarCorr(x))$vcov)
rownames(varcor_random[[i]]) <- c(names(VarCorr(MM_PC_null[[i]][[1]])),
"Residual")
fixef_PC[[i]] <- sapply(MM_PC_null[[i]], fixef)
modmat_fixed[[i]] <- model.matrix(MM_PC_null[[i]][[1]], type = "fixed")
# intercept
XM0 <- modmat_fixed[[i]]
XM0[,-1] <- 0
M0_PC <- XM0%*%fixef_PC[[i]]
M0[[i]] <- M0_PC%*%t(spectra_PCA_loadings)
}
# Effect matrices computation ------------
index <- vector("list", length=length(null_formulas))
fixNames_int <- lapply(Fixvarnames, function(x) c("(Intercept)", x))
colnam <- sapply(modmat_fixed, colnames)
index <- vector("list", length=length(Fixvarnames))
for (i in 1:length(null_formulas)){
id <- c()
for (k in 1:length(fixNames_int[[i]])){
id <- c(id, grep(fixNames_int[[i]][[k]], colnam[[i]]))
}
index[[i]] <- id
}
Mfix <- Mfix_PC <- vector("list", length=length(null_formulas))
names(Mfix) <- fixNames_int
for (i in 1:length(null_formulas)){
XMfix = modmat_fixed[[i]]
XMfix[,-c(index[[i]])] = 0
Mfix_PC[[i]] = XMfix %*% fixef_PC[[i]] # Matrix of the Groupe effect
# backtransform the PC to original coefficients
Mfix[[i]] <- Mfix_PC[[i]]%*%t(spectra_PCA_loadings)
}
######################
# compute the LRT
######################
# objects initialisation ------------
Res_std_error_PC_null <- vector("list", length=length(null_formulas))
loglik_PC_null <- vector("list", length=length(null_formulas))
sumlog <- c()
### compute the LRT ----------------------------
for (i in 1:length(null_formulas)){
Res_std_error_PC_null[[i]] <- sapply(MM_PC_null[[i]], sigma)
# sumlog
loglik_PC_null[[i]] <- sapply(MM_PC_null[[i]], logLik, REML=REML[i])
sumlog[i] <- 2*(sum(loglik_PC_full[i,] - loglik_PC_null[[i]]))
}
names(sumlog) <- names(null_formulas)
sumlog_true <- sumlog
# Graphs ----------------------------
col1="blue"
col2="red"
par(mar=c(4,3,2,6))
dif <- vector(mode = "list")
for (i in 1:length(null_formulas)){
# graphs
mat <- cbind(loglik_PC_null[[i]], loglik_PC_full[i,])
rownames(mat) <- paste0("PC", 1:nPC)
col <- c(col1,col2)
par(xpd=TRUE)
barplot(t(mat), beside=T, ylab="Log-likelihood",
cex.names=0.8, las=2, col=col,
main = paste("log-likelihoods",names(null_formulas)[i]),
xpd=TRUE)
legend("topright", legend = c("loglik restricted","loglik full"),
fill = col, bty = "n", inset=c(-0.1,0))
dif[[i]] <- 2*(mat[,2] - mat[,1])
}
13.3 Booststrap tests
### Bootstrap ---------------------------------------
# bootstrapLT input arguments:
# MM_null = MM_PC_null$volunteer
# useREML=TRUE
# null_formula <- null_formulas[[i]]
pander("running bootstrap ...")
bootstrapLT <- function(useREML, MM_null, null_formula, design, outcomes) {
# simulate y from null models --------
nPC <- length(MM_null)
simulatedY <- c()
for (i in 1:nPC){
ysim <- unlist(simulate(MM_null[[i]], re.form=NA))
simulatedY <- cbind(simulatedY, ysim)
# dim(simulatedY)
}
y <- simulatedY
dimnames(y) <- dimnames(outcomes)
# build restricted model -------
f_null <- parlmer_interaction(design, y, null_formula, REML=useREML)
# build full model -------
f_full <- parlmer_interaction(design, y, form, REML=useREML)
MM_f_null <- f_null$merMod_obj
MM_f_full <- f_full$merMod_obj
# LR -------
loglikelihood_null <- sapply(MM_f_null, logLik, REML=useREML)
loglikelihood_full <- sapply(MM_f_full, logLik, REML=useREML)
ratio <- 2*(loglikelihood_full-loglikelihood_null)
sumlog <- 2*(sum(loglikelihood_full - loglikelihood_null))
return(list(sumlog=sumlog, ratio=ratio))
# returns the summed LLR (sumlog) and the LLR per PC (ratio)
}
# test the function bootstrapLT
bootstrapLT(MM_null = MM_PC_null[["Candies" ]], useREML =REML["Candies" ],
null_formula = null_formulas[["Candies" ]],
design = designInter, outcomes = outcomes)
# Boostrapping: Apply bootstrapLT for each effect
sumlog_boot <- vector("list", length=length(null_formulas))
ratio_boot <- vector("list", length=length(null_formulas))
names(sumlog_boot) <- names(ratio_boot) <- null_effects
set.seed(2018)
for (i in 1:length(null_formulas)){
null_effect <- null_effects[i]
res = replicate(nsim, bootstrapLT(MM_null = MM_PC_null[[null_effect]],
useREML =REML[null_effect],null_formula = null_formulas[[null_effect]],
design = designInter, outcomes = outcomes),
simplify = "array")
sumlog_boot[[i]] <- res["sumlog", ]
sumlog_boot[[i]] <- unlist(sumlog_boot[[i]])
ratio_boot[[i]] <- res["ratio", ]
ratio_boot[[i]] <- do.call(rbind, ratio_boot[[i]])
}
save(sumlog_boot, sumlog_true,ratio_boot, file=file.path(out_path,name_RData))13.4 Boostrap p-values
pval <- c()
for (i in 1:length(null_formulas)){
pval[i] <- (sum(sumlog[i]<sumlog_boot[[i]])+1)/(nsim+1)
}
names(pval) <- names(null_formulas)
pander("p-values")p-values
## Candies Judges CandiesJudges
## 0.0004997501 0.0354822589 0.000499750113.4.1 Graph
difmat <- do.call(cbind, dif)
difmat <- data.frame(PC=substr(rownames(difmat),3,3), difmat)
colnames(difmat) <- c("PC", names(null_formulas))
difmat <- gather(difmat, key=Effect, value = value, Candies , Judges, CandiesJudges)
difmat$Effect <- as.factor(difmat$Effect)
difmat$Effect <- factor(difmat$Effect, levels = c("Candies", "Judges", "CandiesJudges"))
LLRplot <- ggplot(data=difmat, aes(x=PC, y=value, fill = Effect)) +
geom_bar(width=0.5,stat="identity",
position=position_dodge(width=0.5)) +
theme_classic()+labs(title="(Restricted) Log-likelihood Ratios") +
ylab(label="(R)LLR") +
guides(fill=guide_legend(title="Removed effect"))
tab <- paste(c("<", "", "<"), round(pval, 4))
tab <- data.frame(Effect = names(pval), `p-value` = tab,
chi2 = c("<5e-04","-", "-"))
pval_boot.table <- ggtexttable(tab,cols = c("Effect",
"Boostrapped p-value", "Chi2 test"),
rows = NULL)
LLR_pval_plot <- ggarrange(LLRplot, pval_boot.table,
ncol = 1, nrow = 2,
heights = c(1, 0.3), common.legend = TRUE)
# ggexport(LLR_pval_plot, filename = file.path(fig_path,"SDA_LLR_pval_plot.pdf"),
# height = 5, width = 5)
# LLR_pval_plot
# journal article output
p <- LLRplot + annotation_custom(ggplotGrob(pval_boot.table),
xmin = 2, ymin = 100)
p
### Chapter graph
difmat <- do.call(cbind, dif)
difmat <- data.frame(PC=substr(rownames(difmat),3,3), difmat)
colnames(difmat) <- c("PC", names(null_formulas))
difmat <- gather(difmat, key=Effect, value = value, Candies , Judges, CandiesJudges)
difmat$Effect <- as.factor(difmat$Effect)
difmat$Effect <- factor(difmat$Effect,levels(difmat$Effect)[c(2,1,3)])
difmat2 <- cbind(difmat[difmat$Effect=="Candies", "value"],
difmat[difmat$Effect=="Judges", "value"],
difmat[difmat$Effect=="CandiesJudges", "value"])
dimnames(difmat2) <- list(c(1:8), levels(difmat$Effect))
# pdf(file.path(fig_path, "SDA_GLLR.pdf"), width=7,
# height=6, pointsize = 20)
par(mar=c(4,4,4,4), xpd=TRUE)
# plotting settings -------------------------------------------------------
ylim <- range(mat)*c(1,1.5)
angle1 <- rep(c(45,45,135), length.out=7)
angle2 <- rep(c(45,135,135), length.out=7)
density1 <- seq(5,35,length.out=7)
density2 <- seq(5,35,length.out=7)
op <- par(mar=c(4,3,1,1))
barplot(t(difmat2), beside=TRUE,col = gg_color_hue(3),ylab="(R)LLR", xlab="PC",
main = "(Restricted) Log-likelihood Ratios",
angle=angle1[c(2,4,6)], density=density1[c(2,4,6)],
ylim = c(0,200))
barplot(t(difmat2), beside=TRUE, add=TRUE, col = gg_color_hue(3),
ylab="GLLR", xlab="PC",
main = "(Restricted) Log-likelihood Ratios",
angle=angle2[c(2,4,6)], density=density2[c(2,4,6)],
ylim = c(0,200))
legend("topright", legend = c("Candies", "Judges", "Candies*Judges"),
title = "Removed effect:", ncol=1, col = gg_color_hue(3),
fill=gg_color_hue(3),angle=angle1[c(2,4,6)],
density=density1[c(2,4,6)], inset=c(0,0.1), bty="n")
par(bg="transparent")
legend("topright", legend = c("Candies", "Judges", "Candies*Judges"),
title = "Removed effect:",ncol=1, col = gg_color_hue(3),
fill=gg_color_hue(3),angle=angle2[c(2,4,6)],
density=density2[c(2,4,6)],inset=c(0,0.1), bty="n")
# dev.off()
# pdf(file.path(fig_path, "SDA_pval.pdf"))
pval_boot.table
# dev.off()13.5 Histogram drawings
names(sumlog_boot) <- casefold(names(null_formulas), upper = FALSE)
df <- nPC * 4
######################
# Candies
######################
# pdf(file = file.path(fig.path, "hist_time.pdf"), width = 7, height = 5)
par(mar=c(2.1, 2.1, 2.1, 7.5), xpd=TRUE, mfrow=c(1,1))
m=hist(sumlog_boot$candies, freq=F, breaks=100,
xlab="Global Likelihood Ratio Statistic",
xlim=range(sumlog["Candies"], sumlog_boot$candies),
ylim = c(0,0.08),
col = gray67,border = gray67,
main = " Fixed Candies effect", cex.main = 2.2)
lines(density(sumlog_boot$candies), col=darkblue,lwd=2)
lines(dchisq(seq(0,max(sumlog_boot$candies)), df),
col= "limegreen", lwd = 4)
points(sumlog["Candies"], 0, col="red", pch=19, lwd=6)
legend("topright",
legend = c(paste0("True GLLR: ", round(sumlog["Candies"],2)),
"Kernel density",
paste0("chi2 distrib. (df=", df,")")),
col = c("red", darkblue, "limegreen"),
lty=c(NA,1,1),pch=c(19,NA,NA),
inset=c(-0.2,0),box.lty=0, cex = 1.4,
y.intersp = 0.8, lwd=c(4))
# dev.off()
#
######################
# Judges
######################
# pdf(file = file.path(fig.path, "hist_volunteer.pdf"), width = 7, height = 5)
# df*(1/2*dchisq(0:30, 1)+ (1/2*dchisq(0:30, 0)))
par(mar=c(2.1, 2.1, 2.1, 7.5), xpd=TRUE, mfrow=c(2,1))
df <- nPC
m=hist(sumlog_boot$judges, freq=F, breaks=100,
xlab="Global Likelihood Ratio Statistic",
xlim=range(sumlog["Judges"], sumlog_boot$judges),
col = gray67,border = gray67,
main = " Random Judges effect", cex.main = 2.2)
lines(density(sumlog_boot$judges), col=darkblue,lwd=2)
# lines(0:30,(1/2*dchisq(seq(0,30), df)+ (1/2*dchisq(seq(0,30), 0))),
# col= "limegreen", lwd = 4)
points(sumlog["Judges"], 0, col="red", pch=19, lwd=6)
legend("topright",
legend = c(paste0("True GRLLR: ", round(sumlog["Judges"],2)),
"Kernel density"),
col = c("red", darkblue), lty=c(NA,1),pch=c(19,NA),
inset=c(-0.2,0),box.lty=0, cex = 1.4,
y.intersp = 0.8, lwd=c(4))
# dev.off()
######################
# CandiesJudges
######################
# pdf(file = file.path(fig.path, "hist_sampling.pdf"), width = 7, height = 5)
par(mar=c(2.1, 2.1, 2.1, 7.5), xpd=TRUE, mfrow=c(1,1))
m=hist(sumlog_boot$candiesjudges, freq=F, breaks=100,
xlab="Global Likelihood Ratio Statistic",
xlim=range(sumlog["CandiesJudges"], sumlog_boot$candiesjudges),
col = gray67,border = gray67,
main = " Random CandiesJudges effect", cex.main = 2.2)
lines(density(sumlog_boot$candiesjudges), col=darkblue,lwd=2)
points(sumlog["CandiesJudges"], 0, col="red", pch=19, lwd=6)
legend("topright",
legend = c(paste0("True GRLLR: ",
round(sumlog["CandiesJudges"],2)),
"Kernel density"),
col = c("red", darkblue), lty=c(NA,1),pch=c(19,NA),
inset=c(-0.2,0),box.lty=0, cex = 1.4,
y.intersp = 0.8, lwd=c(4))
# dev.off()
# 
names(sumlog_boot) <- casefold(names(null_formulas), upper = FALSE)
df <- nPC * 4
######################
# Candies
######################
# pdf(file = file.path(fig_path, "SDA_hist_Candies.pdf"),
# width = 7, height = 5, pointsize = 15)
par(mar=c(2.1, 2.1, 2.1, 3), xpd=TRUE, mfrow=c(1,1))
m=hist(sumlog_boot$candies, freq=F, breaks=100,
xlab="Global Likelihood Ratio Statistic",
xlim=range(pretty(c(sumlog["Candies"],
sumlog_boot$candies))),
ylim = c(0,0.06),
col = "gray75",border = "gray75",
main = "Candies", cex.main = 2.2)
lines(density(sumlog_boot$candies), col=darkblue,lwd=2,
lty=2)
lines(dchisq(seq(0,max(sumlog_boot$candies)), df),
col= "chartreuse4", lwd = 2)
points(sumlog["Candies"], 0, col="red", pch=19, lwd=6)
legend("topright",
legend = c(paste0("True GLRT: ",
round(sumlog["Candies"],2)),
"Kernel density",
paste0("chi2 distrib. (df=", df,")")),
col = c("red", darkblue, "chartreuse4"),
lty=c(NA,2,1),pch=c(19,NA,NA),
inset=c(-0.1,0),box.lty=0, cex = 1.4,
y.intersp = 0.8, lwd=c(4))
# dev.off()
######################
# Judges
######################
# pdf(file = file.path(fig_path, "SDA_hist_Judges.pdf"),
# width = 7, height = 5, pointsize = 15)
# df*(1/2*dchisq(0:30, 1)+ (1/2*dchisq(0:30, 0)))
par(mar=c(2.1, 2.1, 2.1, 4.5), xpd=TRUE)
df <- nPC
m=hist(sumlog_boot$judges, freq=F, breaks=100,
xlab="Global Likelihood Ratio Statistic",
xlim=range(sumlog["Judges"], sumlog_boot$judges),
col = "gray75",border = "gray75",
main = "Judges",
cex.main = 2.2)
lines(density(sumlog_boot$judges), col=darkblue,lwd=2,
lty=2)
# y <- (1/2*dchisq(seq(0,30), df)+
# (1/2*dchisq(seq(0,30), 0)))
# lines(0:30,y, col= "limegreen", lwd = 4)
points(sumlog["Judges"], 0, col="red", pch=19, lwd=6)
legend("topright",
legend = c(paste0("True GLRT: ", round(sumlog["Judges"],2)),
"Kernel density"),
col = c("red", darkblue), lty=c(NA,2),pch=c(19,NA),
inset=c(-0.2,0),box.lty=0, cex = 1.4,
y.intersp = 0.8, lwd=c(4))
# dev.off()
######################
# CandiesJudges
######################
# pdf(file = file.path(fig_path, "SDA_hist_CandiesJudges.pdf"), width = 7, height = 5, pointsize = 15)
par(mar=c(2.1, 2.1, 2.1, 3), xpd=TRUE, mfrow=c(1,1))
m=hist(sumlog_boot$candiesjudges, freq=F, breaks=100,
xlab="Global Likelihood Ratio Statistic",
xlim=range(pretty(c(sumlog["CandiesJudges"],
sumlog_boot$candiesjudges))),
col = "gray75",border = "gray75",
main = "Candies*Judges", cex.main = 2.2)
lines(density(sumlog_boot$candiesjudges), col=darkblue,lwd=2,
lty=2)
points(sumlog["CandiesJudges"], 0, col="red", pch=19, lwd=6)
legend("topright",
legend = c(paste0("True GLRT: ",
round(sumlog["CandiesJudges"],2)),
"Kernel density"),
col = c("red", darkblue), lty=c(NA,2),pch=c(19,NA),
inset=c(-0.1,0),box.lty=0, cex = 1.4,
y.intersp = 0.8, lwd=c(4))
# dev.off()13.6 Compare the theorerical and true quantiles for fixed effect
13.6.1 Graph histo per PC
# pdf(file = file.path(fig_path, "SDA_hist_perPC_Candies.pdf"),
# width = 11, height = 10, pointsize = 10)
par(mar=c(2, 4, 4, 2), xpd=TRUE, mfrow=c(3,3), cex=1)
ratio_boot_Candies <- ratio_boot$Candies
for (i in 1:dim(ratio_boot_Candies)[2]){
m=hist(ratio_boot_Candies[,i], freq=F, breaks=100,
xlab=paste0("LLR"),
col = "gray75",border = "gray75",
main = paste0("LLR Candies - PC ",i))
lines(density(ratio_boot_Candies[,i]), col=darkblue,lwd=3,
lty=2)
curve(dchisq(x, df=4), col='chartreuse4',
main = "Chi-Square Density Graph",
from=0,to=30, add=TRUE, lwd=3, xpd=F)
# points(sumlog["CandiesJudges"], 0, col="red", pch=19, lwd=6)
}
plot(NULL, xlim=c(0,1), ylim=c(0,1), bty="n", xaxt="n", yaxt="n",
xlab="", ylab="")
legend("topright", legend = c("Kernel density", "chi2 (df=4)"),
col = c(darkblue, "chartreuse4"), lty=c(2,1),
inset=c(0,0),box.lty=0,
y.intersp = 1, lwd=c(4), cex=1.4)
13.6.2 Compare the quantiles by PC and for the GLLR
# By PC
df <- 4
theo_quant <- qchisq(seq(0, 0.9999, 1/nsim), df=df)
# pdf(file = file.path(fig_path,
# paste0("SDA_RLLRqqplot_perPC_","Candies",".pdf")),
# width = 12, height = 6)
par(mfrow=c(2,4), cex=1, mar=c(4,4,4,1))
for (i in 1:dim(ratio_boot_Candies)[2]){
boot_quant <- quantile(ratio_boot_Candies[,i],probs = seq(0, 0.9999, 1/nsim))
qTest <- quantile(ratio_boot_Candies[,i],probs = seq(0.0001,0.999, 1/nsim))
qChi2 <- theo_quant
# boot_quant <- subset(boot_quant, !names(boot_quant) %in% c("0.00%", "99.95%","100.00%"))
plot(boot_quant, qChi2,cex.main=1,
xlim=range(boot_quant),
ylim=range(qChi2),
main = paste("LLR Candies -", colnames(ratio_boot_Candies)[i]), xlab="Sample quantiles",
ylab="theoretical quantiles", cex.main=1)
lines(c(0, min(max(boot_quant), max(qChi2))),
c(0, min(max(boot_quant), max(qChi2))),
col="red", lwd=1.5)
}
# dev.off()
# Globally
df <- nPC*4
Theo_quant <- qchisq(seq(0, 0.9999, 1/nsim), df=df)
boot_quant <- quantile(sumlog_boot$candies,probs = seq(0, 0.9999, 1/nsim))
par(mfrow=c(1,1))
# pdf(file = file.path(fig_path, "SDA_GLLRquantiles_Candies.pdf"),
# width = 6, height = 6, pointsize = 15)
plot(boot_quant,Theo_quant, main = "True and theoretical GLLR for Candies effect", xlab = "Sample quantiles", ylab = "Theoretical quantiles", xlim = range(boot_quant, Theo_quant),
ylim=range(boot_quant, Theo_quant))
lines(x=c(1:max(boot_quant, Theo_quant)), y=c(1:max(boot_quant, Theo_quant)), col="red")
13.6.3 p-value based on chi2 distribution
true GLLR:
## Candies
## 284.9536
p-value
## Candies
## 2.294577e-4213.7 Compare the theorerical and true quantiles for random effects
13.7.1 Graph histo per PC
# pdf(file = file.path(fig_path, "SDA_hist_perPC_Judges.pdf"),
# width = 10, height = 10, pointsize = 10)
par(mar=c(2, 4, 4, 3), mfrow=c(3,3), cex=0.9)
ratio_boot_Judges <- ratio_boot$Judges
for (i in 1:dim(ratio_boot_Judges)[2]){
x <- seq(0,round(max(ratio_boot_Judges[,i])),0.001)
y <- 0.5*(dchisq(x, df=1) + dchisq(x, df=0))
yy <- dchisq(x, df=1)
m=hist(ratio_boot_Judges[,i], freq=F, breaks=100,
xlab=paste0("RLLR"),
col = "gray75",border = "gray75",
main = paste0("RLLR Judges - PC ",i))
lines(density(ratio_boot_Judges[,i]), col=darkblue,lwd=2,
ylim =c(0,7), lty=2)
lines(x,y, type="l", col="red", lwd=2,ylim =c(0,7),xpd=F, lty=3)
lines(x,yy, type="l", col="chartreuse4", lwd=2,ylim = c(0,7),xpd=F)
}
plot(NULL, xlim=c(0,1), ylim=c(0,1), bty="n", xaxt="n", yaxt="n",
xlab="", ylab="")
legend("topright",legend = c("Kernel density",
"mixture of chi2\n(df=0,1)",
"chi2 (df=1)"),
col = c(darkblue, "chartreuse4", "red"), lty=c(2,1,3),
inset=c(0,0),box.lty=0,
y.intersp = 0.8, lwd=c(4),xpd=TRUE)
# dev.off()
par(mfrow=c(1,1))
# pdf(file = file.path(fig_path, "SDA_hist_perPC_CJ.pdf"),
# width = 10, height = 10, pointsize = 10)
par(mar=c(2, 4, 4, 3), mfrow=c(3,3), cex=0.9)
ratio_boot_CJ <- ratio_boot$CandiesJudges
for (i in 1:dim(ratio_boot_CJ)[2]){
x <- seq(0,round(max(ratio_boot_CJ[,i])),0.001)
y <- 0.5*(dchisq(x, df=1) + dchisq(x, df=0))
yy <- dchisq(x, df=1)
m=hist(ratio_boot_CJ[,i], freq=F, breaks=100,
xlab=paste0("RLLR"),
col = "gray75",border = "gray75",
main = paste0("RLLR for CJ - PC ",i))
lines(density(ratio_boot_CJ[,i]), col=darkblue,lwd=2, lty=2)
lines(x,y, type="l", col="red", lwd=2,xpd=F, lty=3)
lines(x,yy, type="l", col="chartreuse4", lwd=2,xpd=F)
}
plot(NULL, xlim=c(0,1), ylim=c(0,1), bty="n", xaxt="n", yaxt="n",
xlab="", ylab="")
legend("topright",legend = c("Kernel density",
"mixture of\nchi2 (df=0,1)",
"chi2 (df=1)"),
col = c(darkblue, "chartreuse4", "red"), lty=c(2,1,3),
inset=c(0,0),box.lty=0,
y.intersp = 0.8, lwd=c(4),xpd=TRUE)
13.7.2 Compare the quantiles
13.7.2.1 Compute theo quant
## [1] 0vect <- c()
for (i in seq(0.0001,0.999,by=1/nsim)){
# print(i)
interval <- c(min(qchisq(i,0),qchisq(i,1)),
max(qchisq(i,0),qchisq(i,1)))
interval <- c(0,10)
vect <- c(vect, uniroot(f,interval,
tol = 0.0001, P = i)$root)
}
theo_quant <- vect
# plot(theo_quant)
### Per PC
############
# bootval <- ratio_boot[["Judges"]]
# bootval <- ratio_boot[["Assessors"]]
# bootval <- ratio_boot[["CandiesJudges"]]
names(ratio_boot) <- c( "Candies","Judges" , "CJ")
for (j in 2:3){# random effects
# pdf(file = file.path(fig_path,
# paste0("SDA_RLLRqqplot_perPC_",names(ratio_boot)[j],".pdf")),
# width = 12, height = 6)
par(mfrow=c(2,4), cex=1, mar=c(4,4,4,1))
for (i in 1:dim(ratio_boot[[j]])[2]){
qTest <- quantile(ratio_boot[[j]][,i],probs = seq(0.0001,0.999, 1/nsim))
qChi2 <- theo_quant
qTest <- subset(qTest, !names(qTest) %in% c("0.00%", "99.95%","100.00%"))
plot(qTest, qChi2,
xlim=range(qTest),
ylim=range(qChi2),
main = paste("RLLR for", names(ratio_boot)[j], "-",
colnames(ratio_boot[[j]])[i]), xlab="Sample quantiles",
ylab="theoretical quantiles", cex.main=1)
lines(c(0, min(max(qTest), max(qChi2))),
c(0, min(max(qTest), max(qChi2))),
col="red", lwd=1.5)
}
# dev.off()
}
14 Session info
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] doBy_4.6.6 car_3.0-8 carData_3.0-4
## [4] mdatools_0.10.3 tidyr_1.1.0 ggpubr_0.3.0
## [7] emdbook_1.3.12 spatstat_1.64-1 rpart_4.1-15
## [10] nlme_3.1-148 spatstat.data_1.4-3 reshape2_1.4.4
## [13] cowplot_1.0.0 ggplot2_3.3.2 gridExtra_2.3
## [16] stringr_1.4.0 pander_0.6.3 plyr_1.8.6
## [19] lme4_1.1-23 Matrix_1.2-18 MBXUCL_0.1
## [22] here_0.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.1.8 splines_4.0.1
## [4] crosstalk_1.1.0.1 digest_0.6.25 foreach_1.5.0
## [7] htmltools_0.5.0 magrittr_1.5 tensor_1.5
## [10] cluster_2.1.0 openxlsx_4.1.5 recipes_0.1.13
## [13] gower_0.2.2 bdsmatrix_1.3-4 colorspace_1.4-1
## [16] haven_2.3.1 xfun_0.15 dplyr_1.0.0
## [19] crayon_1.3.4 jsonlite_1.7.0 survival_3.2-3
## [22] iterators_1.0.12 ape_5.4 glue_1.4.1
## [25] polyclip_1.10-0 gtable_0.3.0 ipred_0.9-9
## [28] R2HTML_2.3.2 webshot_0.5.2 spls_2.2-3
## [31] BiocGenerics_0.34.0 abind_1.4-5 scales_1.1.1
## [34] mvtnorm_1.1-1 rstatix_0.5.0 miniUI_0.1.1.1
## [37] Rcpp_1.0.4.6 xtable_1.8-4 ropls_1.20.0
## [40] foreign_0.8-80 proxy_0.4-24 stats4_4.0.1
## [43] lava_1.6.7 prodlim_2019.11.13 htmlwidgets_1.5.1
## [46] ellipsis_0.3.1 pkgconfig_2.0.3 farver_2.0.3
## [49] nnet_7.3-14 deldir_0.1-25 caret_6.0-86
## [52] tidyselect_1.1.0 labeling_0.3 rlang_0.4.6
## [55] manipulateWidget_0.10.1 later_1.1.0.1 munsell_0.5.0
## [58] phyclust_0.1-29 cellranger_1.1.0 tools_4.0.1
## [61] generics_0.0.2 ade4_1.7-15 pls_2.7-2
## [64] broom_0.5.6 evaluate_0.14 fastmap_1.0.1
## [67] yaml_2.2.1 goftest_1.2-2 ModelMetrics_1.2.2.2
## [70] knitr_1.29 zip_2.0.4 rgl_0.100.54
## [73] purrr_0.3.4 mime_0.9 compiler_4.0.1
## [76] curl_4.3 e1071_1.7-3 ggsignif_0.6.0
## [79] spatstat.utils_1.17-0 tibble_3.0.1 statmod_1.4.34
## [82] stringi_1.4.6 forcats_0.5.0 lattice_0.20-41
## [85] nloptr_1.2.2.1 vctrs_0.3.1 pillar_1.4.4
## [88] lifecycle_0.2.0 data.table_1.12.8 httpuv_1.5.4
## [91] R6_2.4.1 promises_1.1.1 clusterSim_0.48-3
## [94] rio_0.5.16 codetools_0.2-16 boot_1.3-25
## [97] MASS_7.3-51.6 rprojroot_1.3-2 withr_2.2.0
## [100] Deriv_4.0 mgcv_1.8-31 parallel_4.0.1
## [103] hms_0.5.3 timeDate_3043.102 coda_0.19-3
## [106] class_7.3-17 minqa_1.2.4 clValid_0.6-8
## [109] rmarkdown_2.3 bbmle_1.0.23.1 pROC_1.16.2
## [112] numDeriv_2016.8-1.1 Biobase_2.48.0 shiny_1.5.0
## [115] lubridate_1.7.9